Discretization Dependence of Criticality in Model Fluids: a Hard-core Electrolyte 
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Grand canonical simulations at various levels, — 5-20, of fine-lattice discretization are reported 
for the near-critical 1:1 hard-core electrolyte or RPM. With the aid of finite-size scaling analyses it 
is shown convincingly that, contrary to recent suggestions, the universal critical behavior is indepen- 
dent of (^4); thus the continuum ((J" oo) RPM exhibits Ising-type (as against classical, SAW, 
XY, etc.) criticality. A general consideration of lattice discretization provides effective extrapolation 
of the intrinsically erratic ^-dependence, yielding (Tq , pc) ~ (0.04933,0.075) for the = oo RPM. 
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The nature of Coulombic criticality has been an out- 
standing experimental and theoretical issue for more than 
a decade [1,2]. Early experiments suggested that some 
electrolytes exhibit classical or van der Waals (vdW) crit- 
ical behavior while others display Ising-type criticality. 
These results raise the central question: Is Coulombic 
criticality of vdW, Ising, or other type? In contrast to 
simple nonionic fluids, in which experiments, theory, and 
simulation point unequivocally towards Ising behavior, 
the situation in electrolytes is more challenging owing to 
the subtle interplay between strong but screened long- 
range ionic forces and the diverging critical density fluc- 
tuations. 

Nevertheless, recent experiments favor Ising-type crit- 
icality [2,3], as do simulations [4-8] of the simplest 1:1 
equisized hard-core electrolyte model — the so-called re- 
stricted primitive model or RPM [1,2]. On the other 
hand, the recent, most precise simulations [7,8] were per- 
formed on a "finely discretized" or lattice version of the 
model at the relatively low discretization level of C — 5. 
Moreover, in 1998 Valleau and Torrie [9] claimed: "So far 
the results offer no support for the existence of any simple 
Ising-type behavior" in the (more realistic) continuum 
model, while in 2002 they asserted that "the behaviors 
of the continuum and discretized models are strikingly 
dissimilar." The work reported here aims to address this 
specific issue by studying simulations of the RPM at in- 
creasing levels of discretization and, more generally, to 
elucidate discretization effects and to show how one may 
effectively extrapolate to the continuum limit ^ — ^ cxd. 

For completeness, we recall, first, some details of the 
RPM: N equisized hard spheres of diameter a in a volume 
V, half carrying charge -|-go and half —qo, interact via the 
Coulomb potential iq^/Dr in a medium (representing 
solvent) of dielectric constant D. Dimcnsionless reduced 
density and temperature variables are then 



TABLE I. MC estimates of T*{() and p*(C) for the RPM. 



p* = Na^/V = pa^ 



T* = k^TDa/ql 



(1) 



Ref. (C=oo) 


lO^Tc 


10' Pc 


C 


10' Tc 


10' p* 


1996 [4(a)] 


4.87(1) 


6.5(2) 


5 


5.069(2) 


7.90(25) 


1999 [4(b)] 


4.88(2) 


8.0(5) 


8 


4.966(2) 


7.60(20) 


1999 [5(b)] 


4.90(3) 


7.0(5) 


10 


4.952(5) 


7.60(20) 


1999 [6] 


4.92(3) 


6.2(5) 


15 


4.948(5) 


7.55(20) 


2002 [13(c)] 


4.89(3) 


7.6(3) 


20 


4.940(5) 


7.50(20) 


2002 [4(c)] 


4.917(2) 


8.0(5) 


oo 


4.933(5) 


7.50(10) 



The model exhibits phase separation into two neutral 
phases, ion rich and ion poor, at T* 



' ^: see Table I. 



Near criticality, the Debye length — \/T*a? /Airp* 
is small, ~ -ja, and many tightly bound neutral clusters 
form: these cause problems both for theory and simula- 
tion. 

Approximate theories yield classical critical exponents 
[10]; conversely, simulations with finite N and V ex- 
hibit rounded critical points so that finite-size scaling 
techniques are needed to extract reliable conclusions. 
Monte Carlo (MC) simulations through 2001 indicate 

=0.060-0.085: see Table I; however, all these val- 
ues have been derived by assuming Ising-type criticality 
and employing the mixed-field finite-size scaling method 
[11]. Although this approach neglects possibly signifi- 
cant pressure- mixing terms [11(b), 12], the crucial point 
here is that even though these simulations exhibit con- 
sistency with Ising-type criticality, they do not rule out 
other candidates, such as vdW, XY, etc. [7]. 

On the other hand, in recent work, LFP [7(b)] con- 
vincingly resolved Ising-type critical behavior from other 
'nearby' candidates via extensive grand canonical simu- 
lations; but, for computational efficiency, LFP studied 
only the finely discretized, C = 5 version of the RPM, the 
discretization level being defined by ( = a/ao with oq the 
spacing of the 'imposed' lattice [13]. Indeed, the speed of 
such lattice simulations can be 100 times faster than off- 
lattice 'continuum' calculations [13]. Now, when (^oo, 
the discretized model clearly approaches the standard 
continuum RPM; conversely, for C = li it corresponds to 
the most basic lattice model that excludes only double 
occupancy of individual sites. Nevertheless, for C>3, 
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one finds that the discrctizcd models have phase dia- 
grams with gas-liquid separation that are qualitatively 
quite similar to the continuum RPM [13]. 

To conclude that the C = 5 RPM belongs to the (d = 
3)-dimensional Ising (or n = 1) universality class, LFP 
first estimated the critical point precisely using multihis- 
togram reweighting and unbiased finite-size scaling meth- 
ods [14,15]. Then, via unbiased extrapolation techniques 
[14,15], they estimated the critical exponents 7 and v 
obtaining 1.24(3) and 0.63(3). These agree well with 
the Ising values, 7=1.239 and = 0.6303, and exclude 
not only vdW criticality (with 7=1, v=\), but also 
XY (n = 2: with 7 ~ 1.316, v ~ 0.670), self-avoiding walk 
(SAW or n = 0: with 7 1.159, 0.588), and n=l 
criticality with long-range potentials |<p(r)| < $/r^-^. 

Now, if one accepts, say on renormalization-group the- 
oretical grounds, that universal critical behavior is inde- 
pendent of detailed features of a system (as LFP tacitly 
assumed), Ising universality may be considered as estab- 
lished for the continuum RPM. But, in light of [9] specif- 
ically, or more generally, how much faith may be put on 
this presupposition — as yet untested in a complex fluid 
like the RPM? That is the question we answer here, as 
well, incidentally, as obtaining improved estimates of 
and pc for the continuum RPM: see Table I. 

In particular, to support Ising-type criticality in the 
continuum RPM, we have studied the model at dis- 
cretization levels C = 5, 8, 10, 15 and 20 via grand canon- 
ical MC simulations in cubic boxes of dimensions L'^ with 
periodic boundary conditions, and sizes varying from 
L*=L/a = 5 to 12. As seen in LFP, obtaining several 
independent confirmations of the class of criticality for a 
nontrivial system requires a large computational effort. 
On the other hand, the universality class can be deter- 
mined with confidence by evaluating sufficiently precisely 
one universal parameter, say either a critical exponent or 
an amplitude ratio, that distinguishes readily among rea- 
sonable candidates. 

Following LFP [7,15], and previous applications to 
simpler, symmetric systems [16], wc thus focus on the 
grand canonical finite-size parameter Ql defined by the 
dimensionless moment-ratio 

Ql{T; {p)l) = {m:')l/{m^)L. m = p- (p)^, (2) 

where denotes a grand canonical expectation value 
at fixed T and chemical potential adjusted to yield 
the mean density {p)l- As well known, Ql—^\ when 
L — > cx) in any single-phase region, while Ql ^ 1 on the 
coexistence curve diameter, pdiam(T') = \{p^+P~)-: where 
(T) denotes the densities of the coexisting liquid and 
gas phases. 

But the crucial point here is that Ql{Tc,Pc) ap- 
proaches a universal value. Qc, that serves to resolve 
distinct criticality classes rather sharply. Thus for 
Ising systems one has Q]^ = 0.6236 and, as discussed in 



LFP, the classical, SAW and XY values are Qf"^ = 
0.4569 Q^^w = 0, and = 0.8045, while for long- 
range, l/r^^"^ Ising systems, Qc(f) increases almost lin- 
early from vdW to Ising values in the interval | < cr < 
1.966 with Qc(o-=l-9) ^0.600. 

To progress it is necessary to calculate Ql along the Q- 
loci, pq(T; L), defined by the value of {p)l for which Ql 
is maximal at fixed T [7(b), 15]: high precision is essential 
[17]. Then finite-size scaling theory [15] implies that as L 
increases, successive self-intersection points, say T^{L), 
approach the critical point, Tc, rapidly as l/L^^+^^Z'^: see 
Fig. 1. In addition, the difference Ql{T^{L); pq) — Qc 




FIG. 1. Plots of Ql on the Q-loci, pq{T;L), vs. T* for 
C = 5, 8 and 15. The horizontal line marks Ql = Q]?- 

varies, in leading orders, as L~^/'' followed by a jfl'"^''^'^ 
term [15], where J2 is the pressure-mixing coefficient [12]. 

As LFP observed, the self-intersection points T^{L) 
for the C = 5 RPM almost coincide with the universal 
Ising value providing both confirmation of Ising charac- 
ter and a precise estimate of Tc: see Fig. 1. Except for 
a translation by AT* ~ 0.0010, the corresponding plots 
shown in Fig. 1 for C = 8 arc almost identical and pro- 
vide the same results. For 1^ = 15 (as for C, = 1Q and 20, 
not shown) the trends are very similar. Since the val- 
ues Ql'^"^ and are off scale and Qc < 0.600 is im- 
plausible, vdW, XY, and long-range Ising criticality with 
cr < 1.9 are again excluded for all these values of Q. Fur- 
thermore, since the behavior as C, increases by a factor of 
4 changes so little, there seem no grounds to doubt that 
Ising behavior prevails for all C — > 00. 

As a consistency check, consider the slopes, say Qc{L), 
of the plots at the points where Ql = Q]?- by finite- 
size scaling, these should diverge as L^/^. Accordingly, 
in Fig. 2(a) the inverses l/Qc{L) for C==5, 8, 10, 15, 
and 20 are plotted versus L^^^'^ using ^ = tig ~ 0.63 and 
*^vdw = 0.5, along with 1/Q'q{L), the inverse slope at the 
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FIG. 2. (a) Inverse of the slope, 1/Qc, evaluated at 
Ql = Q]^ vs. L~^^* with ■tp = ifig and I'vdW- The solid cir- 
cles arc the inverse of the slopes at the inflection points for 
C = 5. (b) Ratios, Q'c{L;0/Q'c{L;C = 5), vs. L-'^^"'' for 
C = 8, 10, 15 and 20. 



inflection points for C = 5- Evidently, setting ip = VvdW 
does not satisfactorily capture the asymptotic behavior 
of the slopes. On the other hand, the slopes at the in- 
flection points for C = 5 (solid circles) provide definite ev- 
idence for Ising character while the slopes, Q'^, for C>5 
also support Ising criticality. Perhaps the most striking 
fact is that the slopes are so insensitive to (: indeed, 
as seen in Fig. 2(b), the values of Qc{L; () for C>8 are 
no more than 1 or 2% higher than for C = 5. This strik- 
ing independence and the clear verdict of Ising criticality 
for the C = 5 RPM [7,8,15], reinforces the conclusion that 
Ising criticality remains valid in the continuum limit. 

From the plots of Fig. 1, the critical temperature for 
C = 8 can be estimated with the same precision as for 
C = 5. Even though only three system sizes have been 
computed for C=10, 15, and 20, the similar behavior 
seen in Figs. I and 2 for different ^ leads to comparable 
estimates although with larger uncertainties: see Table I. 
As found in [13(c)], T*(C) falls when C increases through 
integral values. The critical densities, p*, can be esti- 
mated by suitably extrapolating the densities pq{Tc] L), 
on the Q-loci at the (estimated) values of T* [7(b), 15]. 
From Table I we see that p*(C) also decreases (in accord 
with [13(c)]). 

Now, to extrapolate effectively to the continuum limit, 
we ask how T*(Q, p^iO^ etc., should vary when ^ in- 
creases [18]. To gain insight consider a d-dimensional 
fluid with pair potential (fi{r) = ^po{r) -t- ^p\{r), where (^o 



is repulsive and of short range, say a, while ipi is smooth, 
attractive and long ranged. If Bi{T) = J dfrfi(r) 
with 1+fi = e-*'i(^)/fcBT^ i = 0, 1, are partial second virial 
coefficients, an approximate, vdW-type equation of state 
is 



p/pkTiT = [1 - pBq{T)]-^ + pBi(T). 



(3) 



On discretization with C,=a/a{), the Bi integrals are re- 
placed by sums, B^ = —^ fi{nao), over integral lat- 
tice vectors, n. One must then ask: How rapidly does B^ 
converge to B°°7 Certainly, the decay of the truncation 
error, £'i(C) = (Bf^ /B^) — 1, when ( oo, will be domi- 
nated by any discontinuities in fi{r) (or its derivatives). 

Specifically, for a hard-core ipo (as needed for the 
RPM) one has B^ /B^ = ¥(0/^(0, where V{C) is the 
volume of a sphere of radius ^ while iV(^) is the number of 
lattice sites satisfying |n| < (J. A heuristic argument [18] 
indicates that Eq should be of rms mag nitude Cd/C(''+^)/^ 
with Cd = 0(1). This is exact for d = 1 and 2 and is sup- 
ported numerically for d = 3 by the scaled plot in Fig. 
3(a) [19]. Evidently, £'o(C) varies wildly and discontinu- 
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FIG. 3. (a) Relative hard-core truncation error, Eo{Q, 
for sc lattices: (a) scaled by C^; (b) magnitude vs. for 
half-odd (x) and integral (o and •) C values. 

ously; and the noisiness persists when ( is restricted to, 
e.g., half-integers: see Fig. 3(b). 

Furthermore, as implied by (3) [18], the erratic behav- 
ior transfers to p*(C) and T*{() and seriously hampers 
extrapolation: see the plots (i) in Fig. 4. We present 
two strategies to mitigate the problem. First, define a 
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FIG. 4. Estimation of pc and for the RPM. Above scale 
breaks: (i) simulation values vs. C,~'^\ below breaks: (ii) vs. 
VC^^(C); (iii) vs. \l(C}—2'f'\ (iv) rescaled simulation values, 
pt(C) andTet(C), vs. C"': see (4). 



modified discretization level C^IO? via Ea[C) = j C} : 
c^ = 25i?o(5) is convenient. Then, as in Fig. 4 plots (ii) 
and, with an e — 2 shift, (iii), examine p*(C), etc., vs. 
Evidently, the behavior is much smoother! 
Second, an effective Bi{T) for the RPM must include 
contributions from the Coulombic interactions. Accord- 
ingly, in the hope of improving convergence, we supple- 
ment Bq by coefficients, bi, independent oi(^, and rescale 
the discretized densities via 



(4) 



and similarly for Tc- Indeed, for the choices b^ — OAB^ 
and bl = 1.7B^ both pJ(C) and T^{C) become almost 
insensitive to see plots (iv). 

From the enhanced plots in Fig. 4 we estimate T* ~ 
0.04933 and 0.075 for the RPM; see Table I. Our 
value for T* agrees well with the (less precise) estimate 
of [6] , but their estimate of p* is very low. Other recent 
estimates of p* encompass our value but the T* estimates 
fall lower [20]. 

In summary, Monte Carlo studies of Ql (T) on the Q- 
loci of the restricted primitive model (RPM) for various 
discretization levels, C = 5-20, have provided convincing 
evidence for Ising-type, as against XY, or SAW, etc., crit- 
icality in the continuum limit. By pinpointing, generally. 



the primary sources of discretization errors [18] we have 
found effective means of estimating precisely the limiting 
critical parameters from data for ( < 10. 
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